In [1]:
import string
import copy
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import os
import sys
import cPickle
import scipy.io as sio
import matplotlib.units
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
from common_functions import *
import half_nanoplate_functions as hnf
import matplotlib_helper_functions as mplhf
%matplotlib inline
#%config InlineBackend.figure_format = 'retina'
#import seaborn as sns
In [2]:
os.chdir("C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Data")
In [3]:
#matplotlib.rcParams['font.family'] = 'Vera'
matplotlib.rcParams['figure.dpi'] = 200.0
#matplotlib.rcParams['axes.labelsize'] = 15
#matplotlib.rcParams['axes.titlesize'] = 18
#matplotlib.rcParams['xtick.labelsize'] = 'large'
#matplotlib.rcParams['ytick.labelsize'] = 'large'
matplotlib.rcParams['mathtext.fontset'] = 'dejavusans'
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['font.sans-serif'] = 'stixsans'
matplotlib.rcParams['legend.handletextpad'] = -0.5
matplotlib.rcParams['legend.numpoints'] = 1
matplotlib.rcParams['legend.scatterpoints'] = 1
matplotlib.rcParams['legend.borderaxespad'] = 0.2
matplotlib.rcParams['legend.fontsize'] = 'small'
matplotlib.rcParams['xtick.labelsize'] = 'small'
matplotlib.rcParams['ytick.labelsize'] = 'small'
matplotlib.rcParams['axes.labelpad'] = 0.5
#matplotlib.rcParams['mathtext.default'] = 'regular'
#matplotlib.rcParams['font.family'] = 'serif'
#matplotlib.rcParams['font.serif'] = 'Times'
#matplotlib.rcParams['text.usetex'] = True
#matplotlib.rcParams['text.latex.unicode'] = True
#matplotlib.rc('text.latex', preamble=r'\usepackage{cmbright}')
matplotlib.rc('text.latex', preamble=[r'\usepackage{DejaVuSans}',
r'\renewcommand*\familydefault{\sfdefault}',
r'\usepackage[T1]{fontenc}',
r'\usepackage{amsmath}'])
In [4]:
# For Matplotlib 2.0
matplotlib.rcParams['xtick.direction'] = 'in'
matplotlib.rcParams['ytick.direction'] = 'in'
matplotlib.rcParams['xtick.top'] = True
matplotlib.rcParams['ytick.right'] = True
In [5]:
store = pd.HDFStore("K:\Pat's_Projects\ParticleTrajectoryData\half_nanoplate_dynamics_processed.h5", mode='r')
keys = store.index['key']
df_l1 = store.get(keys[0])
store.close()
In [20]:
import matplotlib.gridspec as gridspec
um_conv=6.5/60/1.6/2
pmfs_l1 = plt.figure(figsize=(3.34646, 1.8))
df = df_l1.drop_duplicates(['frame', 'track id'])
################
# Plots of PMFs#
################
gs2 = gridspec.GridSpec(1,2)
gs2.update(left=0, right=1.0, bottom=0, top=1)
# pmf leading edge
ax3 = plt.subplot(gs2[0,0])
to_plate_df = df.query('50 < theta < 105')
r_avg = to_plate_df.r.mean()
to_plate_df = to_plate_df.query('@r_avg-10 < r < @r_avg+10')
bin_num = hnf.hist_bin_optimization_continuous(to_plate_df.theta)
hist_to_plate, bins_to_plate = np.histogram(to_plate_df['theta'], bins=bin_num, normed=True)
kernel = scipy.stats.gaussian_kde(to_plate_df['theta'], bw_method=0.05)
theta_range = np.linspace(51,104,400)
pmf = -np.log(kernel.evaluate(theta_range))
plt.plot(theta_range, pmf-np.max(pmf))
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'$-\log(PD)$')
# Annotate Distance of Chords for Optical Binding
indices = scipy.signal.argrelmin(pmf)[0]
angles = theta_range[indices[1:]] - theta_range[indices[:-1]]
chords_um = r_avg*um_conv*2*np.sin(angles*0.5*np.pi/180)
ax3.annotate('', xy=(theta_range[indices[2]],-3.7), xycoords='data',
xytext=(theta_range[indices[3]],-3.7), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='|-|,widthA=0.5,widthB=0.5', shrinkA=0, shrinkB=0))
ax3.annotate('', xy=(theta_range[indices[3]],-3.7), xycoords='data',
xytext=(theta_range[indices[4]],-3.7), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='|-|,widthA=0.5,widthB=0.5', shrinkA=0, shrinkB=0))
ax3.text(theta_range[indices[2]]+(theta_range[indices[3]] - theta_range[indices[2]])*0.2, -3, "{:0.2f}".format(chords_um[2]), va='center', ha='center', fontsize=7)
ax3.text(theta_range[indices[3]]+(theta_range[indices[4]] - theta_range[indices[3]])*0.8, -3, "{:0.2f}".format(chords_um[3]), va='center', ha='center', fontsize=7)
# Add a dotted line where the barrier is
max_indices = scipy.signal.argrelmax(pmf)[0]
plt.plot([theta_range[max_indices[-2]]]*2, [1,-6], 'k--', lw=1.0, zorder=0)
plt.ylim(-5.6,0.25)
# pmf of falling edge
ax4 = plt.subplot(gs2[0,1])
from_plate_df = df.query('240 < theta < 310')
r_avg = from_plate_df.r.mean()
from_plate_df = from_plate_df.query('@r_avg-10 < r < @r_avg+10')
bin_num = hnf.hist_bin_optimization_continuous(from_plate_df.theta)
kernel = scipy.stats.gaussian_kde(from_plate_df['theta'], bw_method=0.05)
theta_range = np.linspace(241,309,400)
pmf = -np.log(kernel.evaluate(theta_range))
plt.plot(theta_range, pmf-np.max(pmf))
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'$-\log(PD)$')
# Annotate Distance of Chords for Optical Binding
indices = scipy.signal.argrelmin(pmf)[0]
angles = theta_range[indices[1:]] - theta_range[indices[:-1]]
chords_um = r_avg*um_conv*2*np.sin(angles*0.5*np.pi/180)
ax4.annotate('', xy=(theta_range[indices[4]],-2.4), xycoords='data',
xytext=(theta_range[indices[5]],-2.4), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='|-|,widthA=0.5,widthB=0.5', shrinkA=0, shrinkB=0))
ax4.text(theta_range[indices[4]]+(theta_range[indices[5]] - theta_range[indices[4]])*0.5, -1.9, "{:0.2f}".format(chords_um[4]), va='center', ha='center', fontsize=7)
# Add a dotted line where the barrier is
max_indices = scipy.signal.argrelmax(pmf)[0]
plt.plot([theta_range[max_indices[3]]]*2, [1,-6], 'k--', lw=1.0, zorder=0)
plt.ylim(-4.1,0.25)
# Label the axes
mplhf.add_axes_label_inches(ax3, (0.068,0.065), '(c)', fontsize=10)
mplhf.add_axes_label_inches(ax4, (0.068,0.065), '(d)', fontsize=10)
# Add labels to each panel
ax3.text(0.5, 0.85, r"at $\pi$/2", transform=ax3.transAxes, ha='center', usetex=True)
ax4.text(0.75, 0.85, r"at 3$\pi$/2", transform=ax4.transAxes, ha='center', usetex=True)
gs2.tight_layout(pmfs_l1, rect=[0, 0, 1, 1], pad=0.05, w_pad=0.1)
#plt.tight_layout(pad=0.01, w_pad=.1)
In [21]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
pmfs_l1.savefig(save_dir+"\Fig_1_pt_2.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
pmfs_l1.savefig(save_dir+"\Fig_1_pt_2.pdf", dpi=800)
Figure 2¶
Electric Field, Calculated Forces, and PMF of Barriers¶
Load the Simulation Intensity and Calculated Potentials¶
This data is from Nishant's simulations and describes the intensity of the electric field over the nanoplate and the potential energy of a particle as it crosses off of the nanoplate. From Nishant:
There are two .mat files. intensity_yz_lcmp.mat contains the field intensities for l=1 to l=5 (fyzsql1 to fyzsql5). yo and zo are the y and z axes in nanometers. The file U_edgesep_lcmp_v2.mat contains the potential barriers. ul0-ul5 are the potentials in units of k_BT and edgesep is the distance from the edge.
In [22]:
efield = mplhf.load_matlab_to_dict("intensity_yz_lcmp_v2.mat")
potential_e = mplhf.load_matlab_to_dict("potentialbarrier_ED_lcmp.mat")
force_e = mplhf.load_matlab_to_dict("potentialbarrier_ED_lcmp.mat")
sample_traj = mplhf.load_matlab_to_dict("sample_traj_l3.mat")
force_along_traj = mplhf.load_matlab_to_dict("Fy_edld_lcmp_fig2.mat")
In [24]:
def running_gaussian_smoothing(x, y, fwhm):
"""This function will smooth the y values with a gaussian weight for
width in the x values
This function will smooth the y values for each x value given. It will
return the same number of points as what goes into the function.
This function takes two arrays where each part is the x or y part of a
2D trajectory. The function will select each point in x and weight all
points in x with a gaussian centered at point x_i. This weighting is
applied to the points in y to do a weighted average of the y values.
:param x: An array of the x values of a trajectory
:param y: An array of the y values of a trajectory. These values will be
averaged with other y values that are nearby in x. The weighting of the y
values will be gaussian with respect to their distance from each x point
:param fwhm: The full-width-half-maximum of the gaussian function in x used
to weight the values in y
"""
new_y = np.zeros(len(x))
std = fwhm/2.3548
for num, x_point in enumerate(x):
gauss = lambda x: np.exp(-(x - x_point)**2/(2 * std**2))
x_weights = gauss(x)
y_weighted = x_weights * y
new_y[num] = sum(y_weighted)/sum(x_weights)
return new_y
In [25]:
from cycler import cycler
import mpl_toolkits.axes_grid.inset_locator
efield_potential = plt.figure(figsize=(3.34646, 5))
# Electric Field Plot
ax1 = plt.subplot(311)
plt.pcolor(efield['ye'], efield['ze'], efield['fyzsql1'].T, cmap='inferno', rasterized=True)
high_int_args_z = np.argmax(efield['fyzsql1'].T[-70:,:110], axis=0)
plt.plot(efield['ye'][:110], efield['ze'][-70+high_int_args_z], 'k', ls='--')
plt.plot(sample_traj['yp'], sample_traj['zp'], 'blue')
plt.plot((-700, 0), (0, 0), '--', color='#FFD700')
plt.plot((0,0), (0, 200), '--', color='#FFD700')
plt.ylabel('z (nm)')
plt.xlabel('Distance from Edge (nm)')
plt.yticks(np.arange(200,-801, -200))
plt.xlim(-652, 590)
# Force along high intensity region
ax2 = plt.subplot(312)
for num in range(0,6):
plt.plot(force_e['edgesep']*10**9, force_e['fy'+str(num)]*10**12, label='$l$='+str(num), ms=2)
plt.plot([-1000,0], [0,0], 'k--', lw=0.5)
plt.xlim(-652, 590)
plt.xlabel('Distance from Edge (nm)')
plt.ylabel('$F_y(l)$ (pN)', usetex=True)
ax2_legend = plt.legend(loc=(0.3, 0.48), labelspacing=-0.1, handletextpad=0.4, handlelength=0.5, frameon=False, borderpad=0.5)
# Make the legend rendered in Latex
ax2_legend_text = ax2_legend.get_texts()
for text_obj in ax2_legend_text:
text_obj.set_usetex(True)
# Inset Force - Force L=0
ax2_inset = mplhf.add_inset_axes_coors(efield_potential, ax2, [0.75, 0.3, 0.35, 0.7])
l_0_min = min(force_e['fy0'])
for num in range(0,6):
plt.plot(force_e['edgesep']*10**9, (force_e['fy'+str(num)]-force_e['fy0'])*10**12)
# Force along trajectory
ax3 = plt.subplot(313)
plt.plot(0,0,label='$l$=0') # Cycle first color
for l in range(1, 6):
y = force_along_traj['fyl'+str(l)]
y = y*10**12
x = force_along_traj['ysep'+str(l)]
new_y = running_gaussian_smoothing(x, y, 12)
plt.plot(x, new_y, label=r'$l$='+str(l))
plt.xlabel('Distance from Edge (nm)')
plt.ylabel('$F_y(l)$ (pN)', usetex=True)
plt.xlim(-652, 590)
plt.plot([-1000,600], [0,0], 'k--', lw=0.5)
mplhf.add_axes_label_inches(ax1, (0.068,0.065), '(a)', fontsize=10, color='white')
mplhf.add_axes_label_inches(ax2, (0.068,0.065), '(b)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.068,0.065), '(c)', corner='upper right', fontsize=10)
plt.tight_layout(pad=0.1, h_pad=0.3, w_pad=0, rect=[0.05, 0.01, 0.99, 0.99])
In [26]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
efield_potential.savefig(save_dir+"\Fig_2.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Manuscript"
efield_potential.savefig(save_dir+"\Fig_2.pdf", dpi=800)
The data for the cumulative dwell time is from:
Ana_16071101 - Hummer Szabo Fit to k(F) vs F.ipynb
And includes the cumulative dwell times as well as the fits to the cumulative dwell times.
The data for the force distributions is from:
Ana_16101701 - Distributions of Stoke's Force from Different Methods
The force versus rate plot was calculated in Ana_16101702 - Hummer Szabo Fit to k(F) vs F using Theta Stokes Force Calc
In [33]:
import cPickle
# From Ana_16071107
f = open("cum_dwell_times.pkl", "r")
cum_dwell_times, fit_dwell_times = cPickle.load(f)
f.close()
# From Ana_16101701
f = open("force_distributions_expt.pkl", 'r')
force_dist = cPickle.load(f)
f.close()
# From Ana_16101702
f = open("bell_hum_szabo_fit.pkl", "r")
bell_hum_szabo_graph = cPickle.load(f)
f.close()
In [35]:
from scipy.optimize import curve_fit
cum_dwell_times_fig = plt.figure(figsize=(3.34646, 4.5))
gs_cum_dwell_force = plt.GridSpec(3, 2, bottom=0.3)
gs_hs_fit = plt.GridSpec(1,1, top=0.2)
# Plot the Cumulative Dwell Times
ax1 = plt.subplot(gs_cum_dwell_force[0,0])
ax3 = plt.subplot(gs_cum_dwell_force[1,0])
ax5 = plt.subplot(gs_cum_dwell_force[2,0])
for num, ax in enumerate([ax1, ax3, ax5], 1):
L = 2*num - 1
ax.plot(cum_dwell_times['dwell_l'+str(L)], cum_dwell_times['cdf_l'+str(L)], 'o', color='#2ca02c')
ax.plot(fit_dwell_times['dwell_fit_l'+str(L)], fit_dwell_times['cdf_fit_l'+str(L)], 'k')
ax.set_ylabel('$1-\mathregular{CDF}$', fontsize=9.5)
ax.set_ylim(0, fit_dwell_times['cdf_fit_l'+str(L)][0]*1.10)
ax5.set_xlabel('Dwell Time (sec)', fontsize=9.5)
ax5.set_ylim(0,1)
ax3.set_xlim(0,0.46)
ax1.set_xticks(np.arange(0, 5, 2))
ax1.set_yticks(np.arange(0, 0.41, 0.2))
ax3.set_xticks(np.arange(0, 0.50, 0.2))
ax3.set_yticks(np.arange(0, 1.1, 0.5))
ax5.set_xticks(np.arange(0, 0.13, 0.06))
ax5.set_yticks(np.arange(0, 0.91, 0.3))
# Plot the Force Distributions
ax2 = plt.subplot(gs_cum_dwell_force[0,1])
ax4 = plt.subplot(gs_cum_dwell_force[1,1])
ax6 = plt.subplot(gs_cum_dwell_force[2,1])
gaussian = lambda x, x0, sigma: (1/(2*np.pi*sigma**2)**0.5)*np.exp((-(x - x0)**2)/(2*sigma**2))
for num, ax in enumerate([ax2, ax4, ax6]):
data = force_dist[num*2]
bin_num = hnf.hist_bin_optimization_continuous(data)
values, bins, patches = ax.hist(data, bins=bin_num, normed=True, color='#2ca02c')
ax.set_ylabel('PDF', fontsize=9.5)
bins_midpoints = (bins[:-1] + bins[1:])/2.0
params, jac = curve_fit(gaussian, bins_midpoints, values, p0=[data.mean(), data.std()])
fit_x = np.linspace(-0.10, 0.5, num=300)
fit_y = gaussian(fit_x, *params)
ax.plot(fit_x, fit_y, 'k')
# print "meand from gauss="+str(params[0])
# print "mean from dist="+str(data.mean())
# print "std="+str(params[1])
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax4.get_xticklabels(), visible=False)
ax2.set_xticks(np.arange(-0.10, 0.5, 0.10))
ax4.set_xticks(np.arange(-0.10, 0.5, 0.10))
ax6.set_xticks(np.arange(-0.10, 0.5, 0.10))
ax2.set_xlim(-0.10, 0.4)
ax4.set_xlim(-0.10, 0.4)
ax6.set_xlim(-0.10, 0.4)
ax6.set_yticks(np.arange(0, 10, 2))
ax2.set_yticks(np.arange(0, 11, 5))
ax4.set_yticks(np.arange(0, 11, 5))
ax6.set_xlabel('Force (pN)', fontsize=9.5)
ax7 = plt.subplot(gs_hs_fit[0,0])
plt.plot(bell_hum_szabo_graph['avg_force'], bell_hum_szabo_graph['k_f'], 'o', color='#2ca02c')
#plt.errorbar(mean_force_list[:-1]*10**12, k_f[:-1], xerr=std_dev_force_list[:-1]*10**12, linestyle='none')
plt.plot(bell_hum_szabo_graph['bell_x'], bell_hum_szabo_graph['bell_y'], 'C3', ls='--', lw=1.2)
plt.plot(bell_hum_szabo_graph['hum_szabo_x'], bell_hum_szabo_graph['hum_szabo_y'], 'k', lw=1.2)
plt.yscale('log')
plt.xlabel('Applied Force (pN)')
plt.ylabel('k(F)')
plt.ylim(10**-1,10**2)
plt.xlim(0, 0.2)
plt.xticks(np.arange(0, 0.25, 0.05))
# bell_eqn_str = r'$k(F) = k_0 e^{Fx^\ddagger}$'
# hum_szabo_eqn_str = r"\begin{align*} k(F) &= k_o\left(1 - \frac{\nu F x^\ddagger}{\Delta G^\ddagger}\right)^{\frac{1}{\nu}-1} \\ &\times e^{\beta \Delta G^\ddagger\left[1 - \left(1 - \frac{\nu F x^\ddagger }{ \Delta G^\ddagger}\right)^{\frac{1}{\nu}}\right]} \end{align*}"
ax7.text(0.01, 1.5*10**0, 'BE', fontdict={'fontsize':'small', 'color':'C3'}, ha='left')
ax7.text(0.01, 2*10**-1, 'DHS', fontdict={'fontsize':'small'}, ha='left')
# ax7.text(0.01, 1*10**1, bell_eqn_str, usetex=True, fontdict={'fontsize':11, 'color':'C3'}, ha='left')
# ax7.text(0.055, 2*10**-1, hum_szabo_eqn_str, usetex=True, fontdict={'fontsize':11}, ha='left')
mplhf.add_axes_label_inches(ax1, (0.080,0.05), '(a)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.080,0.05), '(b)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax5, (0.080,0.05), '(c)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax2, (0.080,0.05), '(d)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax4, (0.080,0.05), '(e)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax6, (0.080,0.05), '(f)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax7, (0.080,0.015), '(g)', corner='upper left', fontsize=10)
# To make more space adjust the font size of the tick labels
for ax in [ax1, ax2, ax3, ax4, ax5, ax6]:
ax.yaxis.set_tick_params(labelsize=7.5)
ax.xaxis.set_tick_params(labelsize=7.5)
#plt.tight_layout(pad=0, h_pad=0.01, w_pad=-0.2, rect=[0,0,1,1])
gs_cum_dwell_force.tight_layout(cum_dwell_times_fig, pad=0, h_pad=0.03, rect=[0,0.4222,1,1])
gs_hs_fit.tight_layout(cum_dwell_times_fig, pad=0, h_pad=0.02, rect=[0,0,1,0.4])
In [36]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
cum_dwell_times_fig.savefig(save_dir+"\Fig_3.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Manuscript"
cum_dwell_times_fig.savefig(save_dir+"\Fig_3.pdf", dpi=800)
In [48]:
import cPickle
f = open("force_distributions_expt.pkl", "r")
force_dist = cPickle.load(f)
f.close()
# From Ana_17030401
f = open("log_pd_experiment.pkl", "r")
log_pd_ls, edge_dist_ls, theta_ls, edge_theta_ls, avg_r_ls = cPickle.load(f)
f.close()
store = pd.HDFStore("K:\Pat's_Projects\ParticleTrajectoryData\half_nanoplate_dynamics_processed.h5", mode='r')
keys = store.index['key']
np_pos = [store.get(v).drop_duplicates(['frame', 'track id']).copy() for v in keys[:-1]]
store.close()
pmf_345 = mplhf.load_matlab_to_dict("pmf_barrier_kernelfit_L345.mat")
umb_samp_pmf = pd.read_csv("emus_pmfs.txt", sep=' ')
In [49]:
from matplotlib.ticker import AutoMinorLocator
from cycler import cycler
um_conv=6.5/60/1.6/2
k_f_vs_force_and_dist = plt.figure(figsize=(5, 3))
# Plot of -log(PD) for Experimental L=1,3,5
###########################################
ax1 = plt.subplot(221)
curve_colors = ['k', '#d62728', '#1f77b4']
for num, v in enumerate(log_pd_ls):
if (num+1)%2 == 0:
continue
plt.plot(edge_dist_ls[num], v, curve_colors[num/2], label='$l$='+str(num+1), lw=1)
plt.xlabel('Distance from Edge ($\mathregular{\mu}$m)', labelpad=0.5)
plt.ylabel(r'$-\log(PD)$')
ax1_legend = plt.legend(loc=(0.60,0.3), labelspacing=-0.1, handletextpad=0.4, handlelength=0.5, frameon=False, borderpad=0.5)
# Make the legend rendered in Latex
ax1_legend_text = ax1_legend.get_texts()
for text_obj in ax1_legend_text:
text_obj.set_usetex(True)
ax1.xaxis.set_minor_locator(AutoMinorLocator(4))
# Set second axis on top for coordinates in theta
# Note: These values in theta are just an estimate because they
# are different for different L's and different radii. The mean values
# of r and position of edge are used to make this axis
ax1a = ax1.twiny()
mean_r_ls = np.array(avg_r_ls).mean()
mean_edge_theta = np.array(edge_theta_ls).mean()
xlim_ax1 = ax1.get_xlim()
theta_lim = (xlim_ax1/(mean_r_ls * um_conv)) * (180.0/np.pi)
ax1a.set_xticks(np.arange(250, 330, 10))
ax1a.set_xlim(theta_lim + mean_edge_theta)
ax1a.set_xlabel(r'$\mathregular{\theta}$ (degrees)', labelpad=5)
markers=['o', '^', 'D', 's', 'x']
# Potential Energy calculated from ED-LD
########################################
ax2 = plt.subplot(223)
ax2.set_prop_cycle(cycler('color', ['#d62728','#9467bd','#8c564b']))
plt.plot(pmf_345['y_values'], pmf_345['pmf3'], label='$l$=3')
plt.plot(pmf_345['y_values'], pmf_345['pmf4'], label='$l$=4')
plt.plot(pmf_345['y_values'], pmf_345['pmf5'], label='$l$=5')
# Find location of optical binding minimum after barrier to align experiment to
l3_edld_opt_bind_min = pmf_345['y_values'][-25:][np.argmin(pmf_345['pmf3'][-25:])]
# Plot L=3 experimental result
l3_expt_opt_bind_min = edge_dist_ls[2][-200:][np.argmin(log_pd_ls[2][-200:])]
shift_edge_dist_l3 = l3_edld_opt_bind_min - l3_expt_opt_bind_min*1000
plt.plot(edge_dist_ls[2]*1000 + shift_edge_dist_l3, log_pd_ls[2] - max(log_pd_ls[2]), '--')
#plt.plot(edge_dist*10**3 + 175, pmf-max(end_pmf), 'b')
plt.xlim(-600, 500)
plt.xlabel('Distance from Edge (nm)')
plt.ylabel(r'$-\log(PD)$')
ax2_legend = plt.legend(loc=(0.45,0.05), labelspacing=-0.1, handletextpad=0.4, handlelength=0.5, frameon=False, borderpad=0.5)
# Set legend text to use Latex
ax2_legend_text = ax2_legend.get_texts()
for text_obj in ax2_legend_text:
text_obj.set_usetex(True)
ax2_xlim = ax2.get_xlim()
# Plot the Umbrella Sampled PMFs
################################
ax3 = plt.subplot(222)
y = umb_samp_pmf.y*10**3 - 937.5
indx_at_600 = np.argmin(abs(y + 600))
for l in range(6):
pmf = umb_samp_pmf[str(l)] / (4.11*10**-21)
plt.plot(y, pmf-pmf[indx_at_600], label=r'$l$='+str(l))
plt.ylabel('pmf ($\mathregular{k_BT}$)')
plt.xlabel('Distance from Edge (nm)')
ax3_legend = plt.legend(loc='lower left', ncol=2, columnspacing=1, labelspacing=-0.1, handletextpad=0.4, handlelength=0.5, frameon=False, borderpad=0.5)
plt.xlim(ax2_xlim)
# Set legend text to use Latex
ax3_legend_text = ax3_legend.get_texts()
for text_obj in ax3_legend_text:
text_obj.set_usetex(True)
# Plot the Umbrella Sampled PMFs with EDLD PMF
##############################################
ax4 = plt.subplot(224)
y = umb_samp_pmf.y*10**3 - 937.5
indx_at_600 = np.argmin(abs(y + 600))
for l in range(3):
pmf = umb_samp_pmf[str(l)] / (4.11*10**-21)
plt.plot(y, pmf-pmf[indx_at_600], label=r'$l$='+str(l))
plt.ylabel('pmf ($\mathregular{k_BT}$)')
plt.xlabel('Distance from Edge (nm)')
ax4_legend = plt.legend(loc="upper left", ncol=3, columnspacing=0.2, labelspacing=-0.1, handletextpad=0.4, handlelength=0.5, frameon=False, borderpad=0.5)
plt.xlim(-600, 200)
ylim_ax4 = np.array([-50, 40])
plt.ylim(ylim_ax4)
plt.yticks(np.arange(-50, 41, 25))
# Set legend text to use Latex
ax4_legend_text = ax4_legend.get_texts()
for text_obj in ax4_legend_text:
text_obj.set_usetex(True)
ax4_2 = ax4.twinx()
ax4_2.set_ylabel(r'$-\log(PD)$')
plt.plot(pmf_345['y_values'], pmf_345['pmf3']-pmf_345['pmf3'][0], '--', color="k", label='$l$=3')
plt.xlim(-600, 200)
scale_factor = 0.05
plt.ylim(ylim_ax4*scale_factor)
mplhf.add_axes_label_inches(ax1, (0.068,0.065), '(a)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax2, (0.068,0.065), '(b)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.068,0.065), '(c)', corner='upper right', fontsize=10)
mplhf.add_axes_label_inches(ax4, (0.068,0.065), '(d)', corner='upper right', fontsize=10)
plt.tight_layout(pad=0, w_pad=0.1, h_pad=0.1)#, rect=[0,0,0.995,0.99])
In [19]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
k_f_vs_force_and_dist.savefig(save_dir+"\Fig_4.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Manuscript"
k_f_vs_force_and_dist.savefig(save_dir+"\Fig_4.pdf", dpi=800)
In [46]:
potential_e = mplhf.load_matlab_to_dict("potentialbarrier_ED_lcmp.mat")
In [47]:
fdtd_pot_e = plt.figure(figsize=(3.34646, 2))
# Potential Energy calculated from FDTD
ax1 = plt.subplot(111)
for num in range(0,6):
plt.plot(potential_e['edgesep']*10**9, potential_e['ul'+str(num)], label='$l$='+str(num), ms=2)
#plt.xlim(np.min(potential_e['edgesep']*10**9), 200)
plt.ylim(-115, 23)
#plt.xticks(np.arange(-500, 251, 250))
plt.yticks(np.arange(0, -125, -25))
plt.ylabel('Potential ($\mathregular{k_B}$T)')
plt.xlabel('Distance from Edge (nm)')
ax1_legend = plt.legend(loc='lower left', labelspacing=-0.1, handletextpad=0.4, handlelength=0.5, frameon=False, borderpad=0.5)
# Set legend text to use Latex
ax1_legend_text = ax1_legend.get_texts()
for text_obj in ax1_legend_text:
text_obj.set_usetex(True)
plt.tight_layout(pad=0.01)
In [48]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
fdtd_pot_e.savefig(save_dir+"\Fig_S1.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Supplemental_Info"
fdtd_pot_e.savefig(save_dir+"\Fig_S1.pdf", dpi=800)
In [55]:
store = pd.HDFStore("K:\Pat's_Projects\ParticleTrajectoryData\half_nanoplate_dynamics_processed.h5", mode='r')
In [56]:
keys = store.index['key']
keys_to_plot = keys[:2]
axes = []
leading_edge_nanoplate = plt.figure(figsize=(3.34646,3))
for i, key in enumerate(keys_to_plot):
L = store.index[store.index['key']==key].L.values[0]
um_conv=6.5/60/1.6/2
df = store.get(key).copy()
df_unique = df.drop_duplicates(['frame', 'track id'])
au_plate_barrier = df_unique.query('50 < theta < 105')
r_avg = au_plate_barrier.r.mean()
au_plate_barrier = au_plate_barrier.query('@r_avg-10 < r < @r_avg+10')
theta_range = np.linspace(50, 105, 400)
theta_range_kde = np.linspace(51, 104, 1000)
ax1 = plt.subplot(len(keys_to_plot), 2, (i*2)+1)
bin_num = hnf.hist_bin_optimization_continuous(au_plate_barrier.theta)
#hist, bins, patches = plt.hist(au_plate_barrier['theta'], bins=bin_num, normed=True)
hist, bins = np.histogram(au_plate_barrier['theta'], bins=bin_num, normed=True)
bins_mid = (bins[1:] + bins[:-1])/2.0
plt.plot(bins_mid, hist, 'o', ms=2)
kernel = scipy.stats.gaussian_kde(au_plate_barrier['theta'], bw_method=0.05)
plt.plot(theta_range_kde, kernel.evaluate(theta_range_kde), lw=1.5)
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'PD')
#plt.title('PMF L='+str(L))
theta_near_barrier = np.linspace(90, 99, 100)
barrier_deg = theta_near_barrier[np.argmin(kernel.evaluate(theta_near_barrier))]
# ax1a = ax1.twiny()
# convert = lambda x: (x-barrier_deg) * (np.pi/180.0) * r_avg * um_conv
# ax1a.set_xlim(ax1.get_xlim())
# ax1a.set_xticks(ax1.get_xticks())
# ax1a.set_xticklabels(["%.1f" % convert(v) for v in ax1.get_xticks()])
# ax1a.set_xlabel('$\mu$m from edge')
ax2 = plt.subplot(len(keys_to_plot), 2, (i*2)+2)
pmf = -np.log(hist)
bins_mid = (bins[:-1] + bins[1:])/2.0
plt.plot(bins_mid, pmf, 'o', ms=2)
plt.plot(theta_range_kde, -np.log(kernel.evaluate(theta_range_kde)), lw=1.5)
#plt.title('PMF L='+str(L))
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'$-\log(PD)$')
axes += [ax1, ax2]
ax1, ax2, ax3, ax4 = axes
mplhf.add_axes_label_inches(ax1, (0.068,0.065), '(a)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax2, (0.068,0.065), '(b)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.068,0.065), '(c)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax4, (0.068,0.065), '(d)', corner='upper left', fontsize=10)
plt.tight_layout(pad=0.1, h_pad=1, w_pad=1)
plt.show()
In [57]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
leading_edge_nanoplate.savefig(save_dir+"\Fig_S2.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Supplemental_Info"
leading_edge_nanoplate.savefig(save_dir+"\Fig_S2.pdf", dpi=800)
In [58]:
keys = store.index['key']
keys_to_plot = keys[:2]
axes = []
falling_edge_nanoplate = plt.figure(figsize=(3.34646,3))
for i,key in enumerate(keys_to_plot):
L = store.index[store.index['key']==key].L.values[0]
df = store.get(key).copy()
df_unique = df.drop_duplicates(['frame', 'track id'])
edge_barrier = df_unique.query('240 < theta < 310')
theta_range = np.linspace(240,310,400)
theta_range_kde = np.linspace(241,309,1000)
r_avg = edge_barrier.r.mean()
um_conv=6.5/60/1.6/2
# Histogram of Positions
ax1 = plt.subplot(len(keys_to_plot), 2, (i*2)+1)
bin_num = hnf.hist_bin_optimization_continuous(edge_barrier.theta)
hist, bins = np.histogram(edge_barrier['theta'], bins=bin_num, normed=True)
bins_mid = (bins[1:] + bins[:-1])/2.0
plt.plot(bins_mid, hist, 'o', ms=2)
kernel = scipy.stats.gaussian_kde(edge_barrier['theta'], bw_method=0.03)
plt.plot(theta_range_kde, kernel.evaluate(theta_range_kde), lw=1.5)
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'PD')
theta_near_barrier = np.linspace(270, 280, 100)
barrier_deg = theta_near_barrier[np.argmin(kernel.evaluate(theta_near_barrier))]
# PMF
ax2 = plt.subplot(len(keys_to_plot), 2, (i*2)+2)
pmf = -np.log(hist)
bins_mid = (bins[:-1] + bins[1:])/2.0
plt.plot(bins_mid, pmf, 'o', ms=2)
plt.plot(theta_range_kde, -np.log(kernel.evaluate(theta_range_kde)), lw=1.5)
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'$-\log(PD)$')
axes += [ax1, ax2]
ax1, ax2, ax3, ax4 = axes
mplhf.add_axes_label_inches(ax1, (0.068,0.065), '(a)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax2, (0.068,0.065), '(b)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.068,0.065), '(c)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax4, (0.068,0.065), '(d)', corner='upper left', fontsize=10)
plt.tight_layout(pad=0.1, h_pad=1, w_pad=1)
plt.show()
In [59]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
falling_edge_nanoplate.savefig(save_dir+"\Fig_S3.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Supplemental_Info"
falling_edge_nanoplate.savefig(save_dir+"\Fig_S3.pdf", dpi=800)
In [60]:
keys = store.index['key']
keys_to_plot = keys[:2]
axes = []
pdf_pmf_nanoplate = plt.figure(figsize=(3.34646,3))
for i,key in enumerate(keys_to_plot):
L = store.index[store.index['key']==key].L.values[0]
df = store.get(key).copy()
df_unique = df.drop_duplicates(['frame', 'track id'])
edge_barrier = df_unique.query('110 < theta < 238')
theta_range = np.linspace(110,238,400)
theta_range_kde = np.linspace(120,230,1000)
r_avg = edge_barrier.r.mean()
um_conv=6.5/60/1.6/2
# Histogram of Positions
ax1 = plt.subplot(len(keys_to_plot), 2, (i*2)+1)
bin_num = hnf.hist_bin_optimization_continuous(edge_barrier.theta)
hist, bins = np.histogram(edge_barrier['theta'], bins=bin_num, normed=True)
bins_mid = (bins[1:] + bins[:-1])/2.0
plt.plot(bins_mid, hist, 'o', ms=2)
kernel = scipy.stats.gaussian_kde(edge_barrier['theta'], bw_method=0.2)
plt.plot(theta_range_kde, kernel.evaluate(theta_range_kde), lw=1.5)
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'PD')
if i==1:
plt.ylim(0,0.012)
#plt.text(0.5, 1.1, 'Particle Positions L='+str(L), fontdict={'fontsize':'large'}, transform=ax1.transAxes, va='center', ha='center')
# PMF
ax2 = plt.subplot(len(keys_to_plot), 2, (i*2)+2)
pmf = -np.log(hist)
bins_mid = (bins[:-1] + bins[1:])/2.0
plt.plot(bins_mid, pmf, 'o', ms=2)
plt.plot(theta_range_kde, -np.log(kernel.evaluate(theta_range_kde)), lw=1.5)
#plt.text(0.5, 1.1, 'PMF L='+str(L), fontdict={'fontsize':'large'}, transform=ax2.transAxes, va='center', ha='center')
plt.xlabel(r'$\theta$ (degrees)')
plt.ylabel(r'$-\log(PD)$')
axes += [ax1, ax2]
ax1, ax2, ax3, ax4 = axes
mplhf.add_axes_label_inches(ax1, (0.068,0.065), '(a)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax2, (0.068,0.065), '(b)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.068,0.065), '(c)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax4, (0.068,0.065), '(d)', corner='upper left', fontsize=10)
plt.tight_layout(pad=0.1, h_pad=1, w_pad=1)
plt.show()
In [61]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
pdf_pmf_nanoplate.savefig(save_dir+"\Fig_S4.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Supplemental_Info"
pdf_pmf_nanoplate.savefig(save_dir+"\Fig_S4.pdf", dpi=800)
In [62]:
store.close()
Schematic of Barrier Crossing¶
In [18]:
def create_schematic_barrier_crossing(ax1, ax2):
'''This function creates a schematic of the barrier crossing with the
relevant parameters listed. The two arguments are the axis objects that
you would like the schematic plotted on.
:param ax1: Matplotlib axis object for plotting the F=0 schematic of
barrier crossing.
:param ax1: Matplotlib axis object for plotting the F>0 schematic of
barrier crossing.
'''
# Create a curve for the equilibrium condition
x = np.linspace(-10,10,1000)
equil_curve = np.poly1d([-5,0,200,0])
# Find the location of the minimum of the well and set it
# to [0,0]
y_min_equil = np.min(equil_curve(x[:len(x)/2]))
x_min_equil = x[np.argmin(equil_curve(x[:len(x)/2]))]
equil_curve_translated = lambda x: equil_curve(x+x_min_equil)+abs(y_min_equil)
# Find the location of the maximum of the barrier
plt.sca(ax1)
y_max_equil = np.max(equil_curve_translated(x[len(x)/2:]))
x_max_equil = x[len(x)/2:][np.argmax(equil_curve_translated(x[len(x)/2:]))]
ax1.plot(x, equil_curve_translated(x))
plt.ylim(-200, 1500)
plt.xlim(-4.5, 11)
plt.ylabel('Potential Energy')
plt.xlabel('Reaction Coordinate')
plt.setp(ax1.get_yticklabels(), visible=False)
plt.setp(ax1.get_xticklabels(), visible=False)
# Create a curve for the equilibrium condition
x = np.linspace(-10,7.8,1000)
f_curve = np.poly1d([-5,0,125,0])
# Find the location of the minimum of the well and set it
# to [0,0]
y_min_f = np.min(f_curve(x[:len(x)/2]))
x_min_f = x[np.argmin(f_curve(x[:len(x)/2]))]
f_curve_translated = lambda x: f_curve(x+x_min_f)+abs(y_min_f)
plt.sca(ax2)
y_max_f = np.max(f_curve_translated(x[len(x)/2:]))
x_max_f = x[len(x)/2:][np.argmax(f_curve_translated(x[len(x)/2:]))]
ax2.plot(x, f_curve_translated(x))
plt.ylim(-200, 1500)
plt.xlim(-4.5, 11)
plt.ylabel('Potential Energy')
plt.xlabel('Reaction Coordinate')
plt.setp(ax2.get_yticklabels(), visible=False)
plt.setp(ax2.get_xticklabels(), visible=False)
# Making the figure look pretty
# Remove spines and ticks
for ax in [ax1,ax2]:
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.yaxis.set_ticks_position('none')
ax.xaxis.set_ticks_position('none')
#ax.axis['xzero'].set_axisline_style('-|>')
ax.annotate('', xy=(1.05,0), xycoords='axes fraction',
xytext=(0,0), textcoords='axes fraction',
arrowprops=dict(facecolor='black', arrowstyle='-|>',
shrinkA=0, shrinkB=0))
ax.annotate('', xy=(0,1.05), xycoords='axes fraction',
xytext=(0,0), textcoords='axes fraction',
arrowprops=dict(facecolor='black', arrowstyle='-|>',
shrinkA=0, shrinkB=0))
ax1.plot([0,x_max_equil+1], [0,0], 'k', lw=1)
ax2.plot([0,x_max_f+1], [0,0], 'k', lw=1)
ax1.annotate('', xy=(x_max_equil,y_max_equil), xycoords='data',
xytext=(x_max_equil,0), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='<|-|>',
shrinkA=0, shrinkB=0))
ax2.annotate('', xy=(x_max_f,y_max_f), xycoords='data',
xytext=(x_max_f,0), textcoords='data',
arrowprops=dict(facecolor='black', arrowstyle='<|-|>',
shrinkA=0, shrinkB=0))
ax1.text(0, -30, r'$0$', va='top', ha='center')
ax1.text(x_max_equil, -30, r'$x^{\ddag}$', va='top', ha='center')
ax1.text(x_max_equil+0.5, 0.3*y_max_equil, r'$\Delta G^{\ddag}$', va='center', ha='left')
k0_arrow_props = dict(facecolor='black', arrowstyle='-|>', connectionstyle='arc3,rad=-0.5', shrinkA=0, shrinkB=0)
ax1.annotate(r'$k_0$', xytext=[0.35*x_max_equil, 0.8*y_max_equil], textcoords='data', xy=[1.01*x_max_equil, 1.1*y_max_equil], xycoords='data', ha='center', va='center', arrowprops=k0_arrow_props)
ax2.text(x_max_f+0.5, 0.3*y_max_f, r'$\Delta U(F)$', ha='left', va='center')
k0_arrow_props = dict(facecolor='black', arrowstyle='-|>', connectionstyle='arc3,rad=-0.5', shrinkA=0, shrinkB=0)
ax2.text(0.35*x_max_f, 1.40*y_max_f, r'$k(F) > k_0$', ha='center', va='center')
ax2.annotate(r'', xytext=[0.35*x_max_f, 0.8*y_max_f], textcoords='data', xy=[1.3*x_max_f, 1.1*y_max_f], xycoords='data', ha='center', va='center', arrowprops=k0_arrow_props)
ax1.text(0.5, 0.9, '$F=0$', transform=ax1.transAxes, ha='center', va='center')
ax2.text(0.5, 0.9, '$F>0$', transform=ax2.transAxes, ha='center', va='center')
# ax1.annotate('', xy=(1.05,0), xycoords='axes fraction', xytext=(0,0), textcoords='axes fraction', arrowprops=dict(facecolor='black', arrowstyle='-|>', shrinkA=0, shrinkB=0))
# ax1.annotate('', xy=(0,1.05), xycoords='axes fraction', xytext=(0,0), textcoords='axes fraction', arrowprops=dict(facecolor='black', arrowstyle='-|>', shrinkA=0, shrinkB=0))
In [31]:
barrier_crossing_schematic = plt.figure(figsize=(3.34646, 2.0))
ax1 = plt.subplot(121)
ax2 = plt.subplot(122, adjustable='box-forced')
create_schematic_barrier_crossing(ax1, ax2)
plt.tight_layout(pad=0.6)
In [32]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
barrier_crossing_schematic.savefig(save_dir+"\Fig_S5.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Manuscript"
barrier_crossing_schematic.savefig(save_dir+"\Fig_S5.pdf", dpi=800)
Comparison of k(F) vs F for Different Experiments¶
In [31]:
# From Ana_16101702
f = open("bell_hum_szabo_fit.pkl", "r")
bell_hum_szabo_graph = cPickle.load(f)
f.close()
# From Ana_17031301
f = open("curtis_bell_hum_szabo_fit.pkl", "r")
curtis_data = cPickle.load(f)
f.close()
In [32]:
curtis_data.keys()
Out[32]:
In [34]:
compare_curtis_data_fig = plt.figure(figsize=(3.34646, 2))
ax = plt.subplot()
# Plot my data
plt.plot(bell_hum_szabo_graph['avg_force'], bell_hum_szabo_graph['k_f'], 'o', color='#2ca02c')
plt.plot(bell_hum_szabo_graph['bell_x'], bell_hum_szabo_graph['bell_y'], 'C3', ls='--', lw=1.2)
plt.plot(bell_hum_szabo_graph['hum_szabo_x'], bell_hum_szabo_graph['hum_szabo_y'], 'k', lw=1.2)
plt.yscale('log')
plt.xlabel('Applied Force (pN)')
plt.ylabel('k(F)')
# plt.ylim(10**-1,10**2)
plt.xlim(0, 0.2)
plt.ylim(1.5*10**-1, 10**3)
plt.xticks(np.arange(0, 0.25, 0.05))
ax.text(0.01, 1.5*10**0, 'BE', fontdict={'fontsize':'small', 'color':'C3'}, ha='left')
ax.text(0.01, 2*10**-1, 'Dudko-Hummer-Szabo', fontdict={'fontsize':'small'}, ha='left')
# Plot Curtis' data
plt.plot(curtis_data['avg_force'], curtis_data['k_f'], '^', color='b')
plt.plot(curtis_data['bell_x'], curtis_data['bell_y'], 'C3', ls='--', lw=1.2)
plt.plot(curtis_data['hum_szabo_fit_x'], curtis_data['hum_szabo_fit_y'], 'k', lw=1.2)
plt.tight_layout(pad=0.1)
In [35]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
compare_curtis_data_fig.savefig(save_dir+"\Fig_S6.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Supplemental_Info"
compare_curtis_data_fig.savefig(save_dir+"\Fig_S6.pdf", dpi=800)
Comparison of Different Definitions of the Edge¶
In [31]:
# From Ana_17031302
f = open("edge_definition.pkl", "r")
compare_edge_data = cPickle.load(f)
f.close()
In [48]:
compare_edge_fig = plt.figure(figsize=(3.34646, 2.5))
ax1 = plt.subplot(221)
plt.imshow(compare_edge_data['median_image'], cmap='gray')
plt.plot(compare_edge_data['edge_x'], compare_edge_data['edge_y_center'], 'r', lw=1)
plt.xlim(275,390)
plt.ylim(250, 100)
plt.xlabel('x (pixels)')
plt.ylabel('y (pixels)')
ax2 = plt.subplot(222)
plt.plot(compare_edge_data['edge_dist_center']*1000, compare_edge_data['logpd'], lw=1)
plt.xlim(-600, 400)
plt.xlabel('Distance from Edge (nm)')
plt.ylabel(r'$-\log(PD)$')
ax3 = plt.subplot(223)
plt.imshow(compare_edge_data['median_image'], cmap='gray')
plt.plot(compare_edge_data['edge_x'], compare_edge_data['edge_y_above'], 'r', lw=1)
plt.xlim(275,390)
plt.ylim(250, 100)
plt.xlabel('x (pixels)')
plt.ylabel('y (pixels)')
ax4 = plt.subplot(224)
plt.plot(compare_edge_data['edge_dist_above']*1000, compare_edge_data['logpd'], lw=1)
plt.xlim(-600, 400)
plt.xlabel('Distance from Edge (nm)')
plt.ylabel(r'$-\log(PD)$')
mplhf.add_axes_label_inches(ax1, (0.068,0.065), '(a)', corner='upper left', fontsize=10, color='white')
mplhf.add_axes_label_inches(ax2, (0.068,0.065), '(b)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.068,0.065), '(c)', corner='upper left', fontsize=10, color='white')
mplhf.add_axes_label_inches(ax4, (0.068,0.065), '(d)', corner='upper left', fontsize=10)
plt.tight_layout(pad=0.1)
plt.show()
In [49]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
compare_edge_fig.savefig(save_dir+"\Fig_S7.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Supplemental_Info"
compare_edge_fig.savefig(save_dir+"\Fig_S7.pdf", dpi=800)
Plot of $\Delta U$ and $x_0$ vs $l$¶
In [5]:
# From Ana_17032901
f = open("graph_del_g_x_o.pkl", "r")
graph_values = cPickle.load(f)
f.close()
# From Ana_16101702
f = open("bell_hum_szabo_model_fit_params.pkl", "r")
model_params = cPickle.load(f)
f.close()
In [45]:
us_param_fig = plt.figure(figsize=(3.34646, 2))
ax1 = plt.subplot()
plt.plot(graph_values['ls'], graph_values['del_gs'], '-o')
plt.xlabel(r'$l$', usetex=True)
plt.ylabel(r'$\Delta G^\ddagger$ or $\Delta U$ ($k_B T$)', usetex=True, fontdict={'color':'C0'})
print model_params['dhs_n0.5_del_g']
plt.plot(0, model_params['dhs_n0.5_del_g'], "^", color="C0", mfc='none')
plt.plot(0, model_params['dhs_n0.66_del_g'],"o", color="C0", mfc='none')
plt.xticks([0,1,2,3])
ax1_2 = ax1.twinx()
ax1_2.plot(graph_values['ls'], graph_values['x_bs'], '-o', color='C3')
plt.ylabel(r'$x_b$ (nm)', usetex=True, fontdict={'color':'C3'}, labelpad=4)
plt.plot(0, model_params['dhs_n0.5_xdd']*10**9, "^", color="C3", mfc='none')
plt.plot(0, model_params['dhs_n0.66_xdd']*10**9,"o", color="C3", mfc='none')
ax1_2.spines['left'].set_color("C0")
ax1_2.spines['right'].set_color("C3")
plt.setp([ax1.get_yticklines(), ax1.get_yticklabels()], color="C0")
plt.setp([ax1_2.get_yticklines(), ax1_2.get_yticklabels()], color="C3")
plt.tight_layout(pad=0.01, rect=[0.01,0,1,1])
plt.show()
In [46]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
us_param_fig.savefig(save_dir+"\Fig_S11.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Supplemental_Info"
us_param_fig.savefig(save_dir+"\Fig_S11.pdf", dpi=800)
Inverse Guassian Figure¶
In [104]:
'''From Curtis' experiment and simulations'''
expt_dwell = mplhf.load_matlab_to_dict('experimental dwelltimes l6.mat')
expt_dwell = expt_dwell[expt_dwell.keys()[0]]
expt_fit = mplhf.load_matlab_to_dict('theoretical data with experimental velocities.mat')
expt_fit = expt_fit[expt_fit.keys()[0]]
theory_barrier = mplhf.load_matlab_to_dict('theoretical mean6 sigma2 d8 k1.mat')
theory_barrier = theory_barrier[theory_barrier.keys()[0]]
theory_no_barrier = mplhf.load_matlab_to_dict('theoretical mean6 sigma2 d8 kNA.mat')
theory_no_barrier = theory_no_barrier[theory_no_barrier.keys()[0]]
In [121]:
exp_fit = lambda x, p1, p2, p3: p1*np.exp(-p2 * (x + p3))
inverse_gauss_fig = plt.figure(figsize=(3.34646, 5))
calc_d_t = lambda array: 1 - np.cumsum(array)/float(sum(array))
inset_coords = [0.58, 0.3, 0.39, 0.6]
ax1 = plt.subplot(311)
plt.plot(expt_dwell[:,1], expt_dwell[:,3], 'o')
plt.plot(expt_fit[:,0], expt_fit[:,1])
plt.xlabel('Dwell Time (sec)')
plt.ylabel('PD')
plt.xlim(0,8)
ax2 = plt.subplot(312)
plt.plot(theory_barrier[:,1], theory_barrier[:,0], 'o', ms=3.7)
x_data = theory_barrier[:,1][8:]
y_data = theory_barrier[:,0][8:]
params, cov = scipy.optimize.curve_fit(exp_fit, x_data, y_data, maxfev=2000)
x_fit = np.linspace(0, x_data[-1], 1000)
y_fit = exp_fit(x_fit, *params)
plt.plot(x_fit, y_fit)
plt.xlabel('Dwell Time (sec)')
plt.ylabel('PD')
plt.xlim(0,8)
plt.ylim(-0.05, 1.3)
ax3 = plt.subplot(313)
plt.plot(theory_no_barrier[:,1], theory_no_barrier[:,0], 'o', ms=3.7)
x_data = theory_no_barrier[:,1][8:]
y_data = theory_no_barrier[:,0][8:]
params, cov = scipy.optimize.curve_fit(exp_fit, x_data, y_data, maxfev=2000)
x_fit = np.linspace(0, x_data[-1], 1000)
y_fit = exp_fit(x_fit, *params)
plt.plot(x_fit, y_fit)
plt.xlabel('Dwell Time (sec)')
plt.ylabel('PD')
plt.xlim(0,8)
plt.ylim(-0.05, 1.3)
plt.tight_layout(pad=0.1, h_pad=0.9)
# Inset CDF
ax1_inset = mplhf.add_inset_axes_coors(inverse_gauss_fig, ax1, inset_coords)
plt.plot(expt_dwell[:,1], expt_dwell[:,2])
plt.plot(expt_fit[:,0], calc_d_t(expt_fit[:,1]))
plt.ylabel('1-CDF')
#plt.xlabel('Dwell Time (sec)')
# Inset CDF
ax2_inset = mplhf.add_inset_axes_coors(inverse_gauss_fig, ax2, [0.58, 0.21, 0.39, 0.6])
plt.plot(theory_barrier[:,1], calc_d_t(theory_barrier[:,0]))
plt.ylabel('1-CDF')
#plt.xlabel('Dwell Time (sec)')
# Inset CDF
ax3_inset = mplhf.add_inset_axes_coors(inverse_gauss_fig, ax3, [0.58, 0.21, 0.39, 0.6])
plt.plot(theory_no_barrier[:,1], calc_d_t(theory_no_barrier[:,0]))
plt.ylabel('1-CDF')
ax2.text(0.35,0.87, 'With Barrier', transform=ax2.transAxes, ha='center', va='center')
ax3.text(0.35,0.87, 'No Barrier', transform=ax3.transAxes, ha='center', va='center')
#plt.xlabel('Dwell Time (sec)')
mplhf.add_axes_label_inches(ax1, (0.068,0.065), '(a)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax2, (0.068,0.065), '(b)', corner='upper left', fontsize=10)
mplhf.add_axes_label_inches(ax3, (0.068,0.065), '(c)', corner='upper left', fontsize=10)
Out[121]:
In [122]:
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Figures"
inverse_gauss_fig.savefig(save_dir+"\Fig_inverse_gauss.png", dpi=600)
save_dir = "C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Manuscript"
inverse_gauss_fig.savefig(save_dir+"\Fig_inverse_gauss.pdf", dpi=800)
1 - CDF of Force¶
In [20]:
cdf_result = []
for num, force in enumerate(force_dist):
bin_num = hnf.hist_bin_optimization_continuous(force)
hist, bins = np.histogram(force, bins=6)
midpoints = (bins[:-1] + bins[1:])/2.0
start = np.where(bins > 0)[0][0]
cdf = np.cumsum(hist)/float(np.sum(hist))
print midpoints
plt.plot(midpoints, (num+1)*hist/(1-cdf), 'o', label='L='+str(num+1))
plt.yscale('log')
plt.legend()
#cdf_result.append(hist/(1-cdf))
